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Abstract — Time synchronization is an important aspect of 
sensor network operation. However, it is well known tliat synchro- 
nization error accumulates over multiple hops. This presents a 
challenge for large-scale, multi-hop sensor networks with a large 
number of nodes distributed over wide areas. In this work, we 
present a protocol that uses spatial averaging to reduce error 
accumulation in large-scale networks. We provide an analysis to 
quantify the synchronization improvement achieved using spatial 
averaging and find that in a basic cooperative network, the 
skew and offset variance decrease approximately as 1/N where 
N is the number of cooperating nodes. For general networks, 
simulation results and a comparison to basic cooperative network 
results are used to illustrate the improvement in synchronization 
performance. 

I. Introduction 
A. Synchronization and the Scalability Problem 

Many synchronization techniques have been proposed for 
synchronizing sensor networks [1], [2], [3], [4], [5]. These 
techniques all rely on nodes exchanging packets with timing 
information. Using the exchanged timing information, each 
node can then estimate clock offset and maybe clock skew. 
However, all of these traditional synchronization techniques 
suffer from an inherent scalability problem — synchronization 
error accumulates over multiple hops. At each hop, nodes will 
estimate synchronization parameters, but the estimates will 
have errors. Therefore, when these erroneous parameters are 
used to communicate timing information to the next hop, errors 
will further increase. 

This accumulation of error over multiple hops poses a 
problem as sensor networks are deployed over larger and larger 
areas. The number of hops required to communicate across the 
network increases and, thus, the synchronization error across 
the network increases as well. One possible way to avoid the 
scalability problem is to use a few nodes with powerful radios 
to limit the number of hops required to communicate timing 
information across the network. However, this technique does 
not address the fundamental scalability problem of errors 
accumulating over multiple hops. 

In this work, we consider the use of high density networks 
to mitigate the scalability problem. Recent developments [6], 
[7], [8] suggest that future networks may have extremely large 
numbers of nodes deployed over wide areas. The question 
we consider is whether or not the density of future networks 
can be used to address scalability issues that plague existing 
techniques. 
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B. Motivation for Cooperation 

In order to reduce the scalability problem, we need to 
find ways to reduce the synchronization error at each hop. 
There are two primary ways to accomplish this. The first is 
to collect more timing information. With more timing data, 
nodes can generally make a better estimate of clock skew 
and clock offset. For example, RBS and FTSP both let nodes 
collect many timing data points before estimating clock skew 
and clock offset. A timing data point provides a node with 
the time at a reference clock at a specific time in its local 
time scale. With more data points, synchronization error will 
decrease. This idea is essentially doing a time average to 
estimate clock skew and clock offset. However, this is not 
necessarily practical since it would significantly increase the 
time to synchronize and the amount of network traffic. 

The second primary approach is to improve the quality of 
the timing data point. For example, TPSN and FTSP use MAC 
layer time stamping techniques that are more accurate. How- 
ever, we believe that there is a fundamentally new technique 
for improving data point quality that has not been considered 
before. This new idea is to use spatial averaging to improve 
data point quality. In a high density network, we have a large 
number of surrounding nodes. Instead of only doing a time 
average to estimate the clock skew and clock offset, perhaps 
we can also do a spatial average to improve these estimates. 

C. Approach to Cooperation 

We assume the network is setup such that one node, called 
node 1, has the reference clock that all other nodes want to 
synchronize to. Node 1 will communicate timing information 
to the nodes in its broadcast domain, the i?2 nodes. The i?2 
nodes will then communicate timing information to the nodes 
that are another hop out, the nodes. This process continues 
until all nodes are synchronized (Fig. 1). 

Each node in i?^, i > 1, will use information from the Ri-i 
nodes to estimate its clock skew and clock offset. With these 
parameters, the node will be able to send a sequence of m 
pulses that are approximately d seconds apart in the reference 
time scale, where d and m are pre-specified by the protocol. 
All nodes in the Ri set will be attempting to send pulses at the 
same time. However, due to synchronization error, the pulses 
will only be occurring at approximately the same time. Thus, 
any node in the Ri+i set of nodes will observe m clusters of 
pulses instead of just m individual pulses. 

Since each pulse in a cluster represents one Ri node's 
attempt to transmit at some appropriate time in the reference 
time scale, taking the sample mean of the pulse arrival times in 
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I network density improves synchronization performance. 
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Fig. 1. For increasing i, the Ri nodes are progressively farther and 
farther away from reference node 1, the iio node. Each node in the 
set Ri receives its timing information from a group of nodes in 

each cluster allows us to average out some of the error made by 
any one node in Ri. The process of having each node in Ri+i 
use the sample mean of each cluster as a timing data point is 
the key to spatial averaging. Since each Ri+i node can use 
timing information from many surrounding neighbors in Ri, 
we call our technique cooperative time synchronization. Using 
these timing data points, and some additional information from 
the Ri set, every node in can estimate its clock skew 

and offset. Thus, this process can then repeat to synchronize 
the Ri+2 nodes. 

The difficulty in studying this problem is that, generally, 
the m clusters observed by any particular node in will 
be different from the m clusters observed by any other node 
in Ri+i. This is because the clusters observed by a node will 
depend on where this node is located relative to the Ri set 
of nodes. Therefore, to study how cooperation can improve 
synchronization, we approach the problem in two steps. 

First, we set up a basic cooperative network (Type I net- 
work) that is a base case for cooperation. The key assumption 
in a Type 1 network is that all nodes in Ri+i are in the 
broadcast domain of all nodes in Ri. Note that this is a 
generalization of the non-cooperative situation where timing 
information is passed from one node to the next. Fig. 2 
compares the basic cooperative network to a non-cooperative 
network. With the Type I network we analytically quantify 
how the variance of the skew and offset estimates grow with 
increasing hop number 
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node 1 node 2 node 3 node 4 

Fig. 2. Top: A basic cooperative network (Type I network) where 
timing information is communicated from the Ri set of nodes to the 
Ri+i set. This Type I deployment assumes many nodes in each Ri 
set. Bottom: Assuming each Ri set has only one node, we have the 
non-cooperative situation. 

Second, we use the theoretical results from the basic coop- 
erative network to understand the behavior of a general Type II 
network where nodes are uniformly distributed over a circular 
area. Simulation results are used to illustrate that increasing 



D. Other Related Work 

The traditional synchronization techniques describe in [1], 
[2], [3], [4], [5] all operate fundamentally on the idea of com- 
municating timing information from one set of nodes to the 
next. One other approach to synchronization that has recently 
received much attention is to apply mathematical models 
of natural phenomena to engineered networks. A model for 
the emergence of synchrony in pulse-coupled oscillators was 
developed in [9] for a fully-connected group of identical 
oscillators. In [10], this convergence to synchrony result was 
extended to networks that were not fully connected. 

The convergence result is clearly desirable for synchro- 
nization in networks and in [11] theoretical and simulation 
results suggested that such a technique could be adapted to 
communication and sensor networks. Experimental validation 
for the ideas of [9] was obtained in [12] where the au- 
thors implemented the Reachback Firefly Algorithm (RFA) on 
TinyOS-based motes. 

The problem with these emergent synchronization results is 
that the fundamental theory assumes all nodes have nearly the 
same firing period. Results from [11] and [12] show that the 
convergence results may hold when nodes have approximately 
the same firing period, but the authors of [12] explain that 
clock skew will degrade synchronization performance. Since 
we are not aware of any results that provide an extension 
to deal with networks of nodes with arbitrary firing periods, 
our work focuses on synchronization algorithms that explicitly 
estimate clock skew. 

E. Contributions and Paper Organization 

In this paper, we propose a protocol for time synchroniza- 
tion that uses spatial averaging to improve synchronization 
performance. In this work we make the following analysis: 

1) Mathematically quantify the synchronization error for 
the Type I basic cooperative network. 

2) Through simulations, show that increasing node density 
can decrease synchronization error in general networks. 

The results show that if each node can hear a large 
number of neighboring nodes, then nodes can cooperatively 
generate signals that are less noisy and allow for better 
synchronization performance over multiple hops. The fact that 
more cooperating nodes yields better performance means that 
there exists a new trade-off between network density and 
synchronization performance where more nodes provide better 
synchronization. Even though it is possible to achieve better 
synchronization performance by introducing nodes with pow- 
erful radios to synchronize a large-scale network, cooperative 
time synchronization is an effective alternative technique to 
reducing synchronization error across the network without 
requiring special nodes. 

The remainder of the paper is organized as follows. In Sec- 
tion II we set up the general network assumptions and present 
the synchronization protocol in Section III. The analysis and 
simulations of the protocol for a basic cooperative network are 
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presented in Section IV while a study of cooperation in general 
networks is carried out in Section V. We make concluding 
remarks in Section VI. 

II. System Model 

A. Clock Model 

The behavior of each node i is governed by a clock Ci 
that counts up from 0. The introduction of Ci is important 
since it provides a consistent timescale for node i. This is the 
node's local time scale and in synchronization the node tries to 
estimate how its local clock is related to the reference clock. 

We assume that node 1 contains the reference clock and 
every node in the network is to be synchronized to this clock. 
The clock of node 1, ci, will be defined as ci{t) = t where 
t S [0, oo). We now define the clock of any other arbitrary 
node i, Ci, as 

c,{t)^a,{t-Ai) + -^,{t), (1) 

where 

• Ki is an unknown offset between the start times of Ci 
and ci. 

• Qfi > is an unknown constant for each i. 

• is a stochastic process modelling random timing 
jitter. i{t) is a zero mean Gaussian process with inde- 
pendent and identically distributed Gaussian samples with 
mean zero and variance a^, i.e. 7V(0, ct^). We assume 

< oo and note that is defined in terms of the clock 
of node i. 

Note that this linear relationship is valid for short periods of 
time since we do not explicitly model clock drift. 

B. Transmission Model 

Each node i in the network can transmit short pulses p{t) 
for time synchronization. These are short duration pulses, i.e. 
ultra wideband pulses, and for our purposes we consider them 
to be delta functions 5{t). The particular choice of pit) is 
not important. For the purposes of studying cooperative time 
synchronization, we assume a node receiving the pulse can 
uniquely determine a pulse arrival time, pulses sent from 
different nodes do not overlap, and a node seeing multiple 
pulses can identify the different pulse arrival times. Note 
that only minor modifications of the protocol are needed to 
accommodate other types of pulse shapes [14]. 

We assume that each node has a transmission range of R. 
This means that a node j must be within a distance R from 
a transmitting node i in order to hear pulses from node i. 
Note that the assumption of a circular transmission region is 
made only to simplify the illustration of spatial averaging. 
The synchronization protocol proposed in Section III does 
not require this assumption and most of the results in this 
work will hold under more realistic conditions [15], [16]. 
Since we are dealing with sensor nodes who have short 
transmission distances, we further assume that propagation 
delay is negligible. We make this assumption since from [5] 
we know that propagation delay is less than l^s for distances 
up to 300 meters. Some results on cooperation and propagation 
delay are available in [13]. 



III. Synchronization Protocol 

To start synchronization, the reference node, node 1, will 
send a sequence of m pulses that are d seconds apart. Since 
we assume the nodes have impulse radio transmitters, each 
pulse is extremely narrow in time. The values of d and m 
are parameters of the protocol that are established before 
deploying the network so the values are known by all nodes in 
the network. Therefore, in the time scale of node 1 the pulses 
are transmitted at times tq, tq + d, . . . , tq + d{rn — 1), where 
tq is the time at which the synchronization process started. 
Let node 1 be the only element of the Rq nodes. 

Pulse clusters observed by node 21 

/ 
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Tg+dm TQ+d(m+l) 'UQ+d(m+2) 'Ug+d(2m-l) 

Fig. 3. A node 2i in R2 has clock C2i. This node will see clusters 
of pulse arrivals that are transmitted from a group of nodes in the 
set Ri. These clusters arrive at node 2i around times ro + dm, tq + 
d(m + 1), . . . , To + d{2m — 1) in the time scale of node 1, ci. 

The nodes that are in the broadcast domain of node 1 will 
hear this sequence of m pulses. We call these nodes the Ri 
nodes and each node i £ i > 1, will be denoted by node 
li. The vector of pulse arrival times observed by node li will 
be denoted Yi;. Each node li will be able to estimate its clock 
skew since it knows that node 1 transmitted these pulses d 
seconds apart. Each node li will also predict, in its own time 
scale, when times tq + dm, tq + d{m + 1), . . . , tq + d{2m — 1) 
will occur in the time scale of node 1 and transmit m pulses, 
one at each predicted time. This means that each node li will 
transmit a pulse approximately at times tq + dm, tq + d{m + 
1), . . . ,To + d{2m — 1) in the time scale of node 1. When 
the Ri nodes each transmit their sequence of m pulses, the 
nodes that can hear a subset of the Ri nodes, the R2 nodes, 
will observe clusters of pulses around the times tq +dm, tq + 
d{m + 1), . . . , To + d{2m — 1) since each node 2i can hear 
many i?i nodes (Fig. 3). In fact, we make sure each node 2i 
can hear a cluster by requiring the node to observe at least N 
pulses in each cluster If a node 2i sees less than N pulses in 
a cluster, then it will not make observations. Each node 2i, a 
node i S R2, will note the aiTival time of each pulse in the fcth 
cluster, fc = 1, . . . , 771, and take the sample mean of these times 
to be its fcth observation. Node 2i's vector of observations will 
be denoted as Y2i- Using these m observations, any node 2i 
will be able to estimate its clock skew since it knows that 
these observations should be occurring d seconds apart. As 
well, it will be able to predict in its local time scale when 
times To + d{2m),To + d{2m + 1), . . . , tq + d{3m — 1) will 
occur in the time scale of the reference time. Node 2i will 
then transmit a pulse at each of those predicted times. This 
processes will repeat until all nodes in the network have an 
estimate of their clock skew. Notice that the i?i nodes are not 
required to observe N pulses in each cluster since they will 
always only receive a sequence of m pulses from node 1. Node 
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1 can simply broadcast a special packet to its surrounding 
nodes to identify the i?i nodes. An illustration of the process 
can be found in Fig. 1 and note that nodes will remain silent for 
the remainder of the synchronization process after transmitting 
their m pulses. The cooperation occurs when a node ki in Rk, 
k > 1, can take a sample mean of a cluster of pulse arrivals. 

To obtain the clock offset, the Rk nodes will broadcast a 
packet of information to the i?A:-|-i nodes, fc > 0. This packet 
will contain the value of tq and a number q denoting the 
number of hops out from node 1. For example, node 1 will 
transmit the value of tq and g = to the Ri nodes. The Ri 
nodes will then send tq and q = 1 to the i?2 nodes. In general, 
the Rk nodes will send tq and q ~ k to the Rk+i nodes. Any 
node ki will then know that its first observation approximately 
occurred at time to + drjiq in the time scale of the reference 
time, where the value of q is the one received from set Rk-i- 

We now describe how any node ki can estimate its clock 
skew, clock offset, and its m pulse transmission times. From 
(1), we know that there is a linear relationship between the 
reference clock ci and the clock of node ki, Cki- Node ki will 
have a set of m observations denoted by the to x 1 vector 
Yki, where the elements of the vector are ordered from the 
earliest observation time to the latest observation time. Node 
ki will estimate its clock skew as 



and clock offset as 

Afe, = C(H^H)-iH^Yfe, - (to + dm{k - 1)), 
where C = [0 1], C = [1 0] and 
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(4) 



Note that in the calculation of the clock offset Aki, the term 
To + dm{k — 1) is the time in the time scale of ci that node 
ki should receive its first pulse. Node ki has used the To and 
g = fc — 1 parameters sent to it from the Rk-i nodes. Node ki 
will also estimate its own to pulse transmission times using 



Xk(0 -Q(H^H)-1H^Yh, 



(5) 



where = [1 d{m + I)], for / = 0, 1, . . . , m - 1. Xk^{l) 
is the transmission time of node ki's {I + l)th pulse. A 
pseudo-code description is given in Table I. Note that the 
protocol described above is a completely new approach to the 
asymptotic spatial averaging ideas we studied in [17]. 

IV. Type I: Basic Cooperative Networks 

A. Network Setup 

The most basic and fundamental deployment of nodes that 
effectively employs cooperative time synchronization is the 
case where all N nodes at any given hop contribute to the 
signals observed at the next hop. This Type I deployment is 
illusti-ated in the top of Fig. 2 where each set of nodes Ri, 
i > 1, have N nodes. We see that every node in i?j is in the 
broadcast domain of every node in Ri^i. 



TABLE I 

The synchronization protocol for each node ki, k > 1. 



Cooperative Time Sync 

wait for pulse arrivals, at least A'^ 

per cluster; 
while (number of arrival clusters < m) { 

record arrival time of all pulses; 

listen for packet with tq and q values; 

}; 

for each (pulse arrival cluster j) { 
Yki[j] = sample mean of cluster; 

}; 

skew ak, = C{H^H)-'H^Yk,; 

offset Afe, = C(H^H)-iH^Yfc, - (ro + dmq) 

for (1 from to m— 1) { 

transmission time = C;(H^H)"^H^Yfei 

transmit pulse at Xki{l); 

}; 

while (transmitting pulses) { 

send a packet with values ro and q + 1; 

}; 



B. Analysis 

Due to the scalability problem, we would expect synchro- 
nization error to grow as timing information from node 1 (the 
Ro node) is communicated to a node in the Rk set of nodes, 
k > 0. Therefore, it is of particular interest to quantify how 
the variance of the skew and offset estimates change as a 
function of the hop number k. Looking at the structure of 
the basic cooperative network in Fig. 2, we notice that the 
skew and offset estimates at a node ki must be dependent 
only on the estimates made by the nodes in Rk-i since all 
the information at Rk comes from the Rk-i set of nodes. 
Therefore, to understand the synchronization error growth over 
multiple hops, we need the distribution of the estimates made 
by the nodes in Rk as a function of the distribution of the 
estimates made by the nodes in Rk-i- Theorem 1, below, 
provides us with this characterization. 

In the statement of the theorem we use e; to be the column 
vector of all zeros except for a one in the /th position and 
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where H is from (4). Also, aki and Aki are the clock model 
parameters from (1) for node ki. 

Theorem 1: Assume a Type I basic cooperative network. 

(1) Given the distribution of the 2N x 1 vector of estimates 
made by the Rk-i nodes. 



• N{jik-i,^k-i), 



the distribution of the 2N x 1 vector of estimates made by the 
Rk nodes, 

Ifc ~A/'(/2fc,Sfe), 
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is found as follows: 9^ has mean vector 

= S(Ajfc_i+Bfc) 

afci(To + (fc - l)md~ Afci) 

afcw(To + (fc - l)md - Afcjv) 



(7) 



and covariance matrix 



for 



tk = Cov(4) = Sm, + AfcSfc„iAf 



S,„, = (H^H)-iH^Sw,((H^H)-iH^)^ 



(8) 



where 



dm 



1 



a(fc-l)l 



Dfc = ^ 



Bfc ~ Dl; 






a.ki 
ckfci 



ttfe^ 




^(fc-l)Ar 





— i — 



aki 

Qffci 



ttfc^ 




^kN^kN 



The initial conditions are 



Ml 



and El = (T2(jjrH)-i 



aii(To - All) 
an 

ai^(ro - Ai^) 



(2) The skew estimate and offset estimate for node ki can 
be found as 

an = e^(,_i)+24 (9) 



Afci = e^(i„i)+i6'fc - (to + dm{k ~ 1)) (10) 

for i 1,2, . . . ,iV and fc > 1. A 

The proof of Theorem 1 is found in the appendix. Since 
the distribution of Ok is available, the distribution of aki and 
Afci can be found. In fact, the variance of aki can be found in 
element (2(i-l) + 2, 2(i-l)+2) of Efe in (8) and the variance 
of Afci can be found in element (2(i — 1) + 1, 2{i — 1) + 1). 
The mean of aki is the {2{i — 1) + 2)th element of /Ifc in (7) 
and the mean of Afci can be found from the (2(i — 1) + l)th 
element shifted by tq + dra{k — 1). 

From the statement of Theorem 1, we see that the dis- 
tribution of the estimates made by the Rk nodes, 6k, is 
completely determined from the distribution of ^fc_i. This 
recursive nature comes from the fact that the parameters 
estimated by the Rk nodes is only dependent on the estimates 
made by the i?fc_i nodes. The relationship between ^fc_i and 
Ok can be intuitively understood in two steps. First, ^fc_i 
is the vector of synchronization parameters estimated by the 
nodes in Rk-i- Therefore, these estimates will establish the 
synchronization parameters for the Rk nodes since the i?fc_i 
nodes communicate timing information to the Rk nodes. The 
synchronization parameters for Rk are found as 



AkO 



feCfe_i 



B 



k- 



Second, the Rk nodes will use the timing information from the 
i?fc_i nodes to make an unbiased estimate of the parameters 
dk, which gives us 9k- 

Since any node fcz's skew and offset estimates are found as 
affine transforms of dk in (9) and (10), respectively, we see 
that any estimation errors made by the i?fc_i nodes will be 
propagated to the Rk nodes' estimates of clock skew and clock 
offset. However, the intuitive understanding of cooperative 
time synchronization comes from realizing that the matrix Afc 
takes an "average" over 9k-i thus mitigating the errors made 
by any particular node {k—l)i. As a result, the synchronization 
parameters communicated to the Rk nodes will be less noisy 
and, therefore, the skew and offset estimates made by a node 
ki will have less error. We would, thus, expect the variance of 
the estimates to decrease with increasing N. Notice that our 
Type 1 network analysis does not explicitly utilize the circular 
transmission region with radius R. 

C. Simulation Results 

In Fig. 4 we illustrate the MATLAB simulation results for 
two 20 hop networks, one with N = 2 and the other with 
N = 4. The following parameters were used: 



R = 1 



0.01 



For each network, a set of iV = 20iV + 1 nodes were 
first placed in a Type 1 network deployment. Each node's 
skew parameter was then generated using ai = \Xi\ for 



6 



Xi ~ A/^(l, 0.005), independently for each node i. Node 1 
was assumed to have ai ~ 1. The cooperative time synchro- 
nization protocol was then run 5000 times using the deployed 
network. At each hop, the 5000 skew and offset estimates of 
one chosen node were used to generate the simulated skew and 
offset estimate variance curves shown in Fig. 4. The theoretical 
variance value of the chosen node at each hop was computed 
using the recursive expression found in (8). 



Clock Skew Estimate 




Fig. 4. Var(afci) is plotted in the top figure and Var(Afei) is plotted 
in the bottom figure as a function of k. 

In Fig. 4, we first clearly see that the simulated skew and 
offset variance values nicely match the predicted theoretical 
variance values. As well, the expected decrease in skew and 
offset variance as N increases from 2 to 4 is immediately 
noticeable. In fact, in both the skew variance and offset vari- 
ance curves, we have an approximate halving of the variance 
values as we double N from 2 to 4. Also expected, is that the 
variance values at each hop depend on the particular values 
of a;, z = 1, . . . , A^. This dependence on the ai values result 
in the jagged skew and offset variance curves seen in Fig. 4. 
The N = 2 network had a, values ranging from 0.9073 to 
1.1342, while the TV = 4 network had skew values ranging 
from 0.8339 to 1.1669. 



The problem with having the variance curves depend on 
the actual skew values is that the exact performance of 
cooperative time synchronization is dependent on the network 
realization. However, we find that for ai values that are close 
to and centered around 1, the variance curves follow the trend 
established by the theoretical variance curves for ai — 1, all 
i. This can be seen in Fig. 4 where we have also plotted 
the theoretical curves using ai = 1 for all i for N = 2 
and TV = 4. As a result, the situation where ai = 1, all i, 
can be used to study the the performance improvement of 
cooperative time synchronization without dealing specifically 
with the skew values of individual nodes. 

Therefore, to get a better understanding of how cooperative 
time synchronization improves synchronization performance, 
let us simplify the recursive expression in (8) for the special 
case where = 1 for all i and find a non-recursive expression 
for skew and offset variance. The first thing to note is that 
under the assumption of ai = 1 for all i, Ak ~ A and = 
Yirii are no longer dependent on k. Therefore, writing out the 
recursive expression for (8), we have 



k-2 



(11) 



4=0 



Using (11), Corollary 1 gives us the non-recursive expression 
for skew and offset variance. 

Corollary 1: For a basic cooperative network with a; = 1, 
all i, ctki and Aki have the following mean and variance; 

E{ak^) = 1 



ki 



Var(afei) 

Var(A;,,) 
+ (fc 



12(7^ 



(P{'m — l)'m{m + 1) 



2(fc-l) 
N 



(12) 



2cr2(2m- 1) 

m(TO + 1) 

n9 / 12 
1)' 



(to 

2){k- l)(2fc 



1 



a 

'n 



3) 



2 r 



4(fc- 1)(2to, - 1) 



to(to 
12to 



1 



{m - l){m + 1) 

12to 
{m - 1)(to + 1) 



(13) 



where fc is a positive integer A 

The proof of Corollary 1 is omitted since it is a direct 
simplification of (11). Note that the skew and offset variance 
expressions are only a function of k and not i. The theoretical 
skew and offset variance of the ith node at the fcth hop (node 
ki) can be found in elements (2(i — 1) + 2, 2(z — 1) + 2) and 
(2(i-l) + l, 2(z-l) + l), respectively, oftk in (11). However, 
the skew variance values in elements (2(i — 1) + 2, 2(i — 1)+2), 
i = 1, . . . , iV, are all equal and the offset variance values in 
elements {2{i - 1) + 1, 2(j - 1) + 1), i = 1, . . . , iV, are also 
equal when we assume that = 1 for all i. As a result, we 
can consider the skew and offset variance at a hop k without 
specifying a particular node. Notice also that, besides the sign 
change in the mean of the offset estimate, the skew and offset 
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estimates are unbiased estimates of the clock parameters of 
node ki. 

Looking at the skew and offset variance curves in (12) and 
(13), respectively, we see that the variance growth decreases 
Uke 1/iV. This 1/N factor in both (12) and (13) is expected 
since every node takes the sample mean of TV pulses to be 
an observation. The variance of the observation decreases like 
1 /N because it is a sample mean and, thus, it is not surprising 
that the skew and offset variance values also approximately 
decrease like 1/iV. 

V. Type II: General Networks 
A. Network Setup 

Nodes will not generally be clustered together as in a 
basic cooperative network, but be deployed in a more random 
manner. As a result, to study general network deployments, we 
will consider a Type II situation where nodes are uniformly 
deployed with density p over a circular region of radius LR 
with node 1 at the center. In such a setup, at any hop k, k > 2, 
a node ki in the Rk nodes will see at least N nodes from the 
Rk-i set of nodes. However, the exact number of observed 
nodes will depend on node fcz's location in the region occupied 
by the Rk nodes. 

An illustration of a Type II deployment is shown in Fig. 5. 
We note that the Rq node (node 1) is placed at the center of 
the disk and the Ri nodes occupy a circular region of radius 
R. However, the region occupied by the Rk nodes for k > 2 
is a ring centered around node 1 with a ring thickness of 
dmax,k- For increasing k, the distance from node 1 to the 
inner circular boundary of the region occupied by the Rk set 
of nodes increases. 



network. The primary change that would occur in the analysis 
is the determination of the affine transform 



1 — Afe+i^fc + Bfc 



fc+i- 




Fig. 5. A Type II network deployment. Nodes are deployed with 
uniform density p and node 1 is at the center of the network. 



B. Analysis 

To study a Type II network, we could carry out an analysis 
similar to the one we did for the Type I basic cooperative 
network. Assuming we know the location of all nodes for a 
given network deployment over the circular region of radius 
LR, we would be able to determine the neighbors of each 
node and then readily extend the Type I analysis to this Type II 



However, there are two issues that arise in determining the 
transform matrix A^+i and vector Bfe+i. 

First, since Rk and Rk+i will most likely have different 
numbers of nodes, we immediately see that A^+i will be a 
2|i?fc+i I X 2|i?fc| matrix and B^+i will be a 2|_Rfe+i | x 1 vector, 
where \Rk\ is the cardinality of set Rk- This means that the 
length of vector Ok will change with every hop. 

Second, for any node {k+l)i in Rk+i, the set of cooperating 
nodes in Rk will be different. Thus, A^. will also reflect this 
difference. Therefore, every time we move from hop fc to fc+1, 
the correlation structure of 9k will change. 

Together, these two points suggest that even though it is pos- 
sible to carry out the full analysis, the complexity would make 
the resulting expressions depend on the particular network 
realization and not provide significant insight into the problem. 
In fact, it would be nearly impossible to visualize the result 
without carrying out a numerical evaluation. Since our goal 
is to comprehend the impact of spatial averaging on general 
networks, we choose to proceed directly with simulations and 
compare the results with our analytical expressions for Type 
I networks. 

In the following analysis, we develop a basic understanding 
of what we would expect to see in the simulation results that 
are presented in Section V-C. We assume that the number of 
nodes in any given area of the Type II network is proportional 
to the area. The reason is that for uniformly deployed nodes 
with density p, the average number of nodes in an area A is 
Ap. Note that even though the analysis and simulation results 
for Type II networks use the assumption of a circular trans- 
mission range of R, the simulation results in Section V-C still 
provide valid insight when realistic transmission regions [15], 
[16] are assumed since the figures illustrate synchronization 
error as a function of hop number Therefore, regardless 
of the shape of the transmission region, a node at hop k 
will have received the appropriate synchronization information 
and, thus, our simulation results reflect its synchronization 
performance. 

1 ) Estimation of L: Our first consideration is to estimate 
the number of hops, L, required to communicate timing 
information from node 1 to the edge of the network a distance 
LR away. In order to do this, we need a way to quantify 
dmax,k- In Fig. 6, we illustrate dmax,2 and see that d,nax.2 
is determined by having the intersection of the two radius R 
circles contain an average of N nodes. This is because if we 
increase d„iax,2, then nodes at this increased distance will not 
see N nodes on average and, thus, not be considered an i?2 
node. However, dmax,k > d„iax,2, for k > 2, because the ring 
occupied by the Rk nodes increases in size for increasing k. As 
a result, we choose to be conservative and let d„iax = drnax,2 
approximate dmax,k for all k. This means that our estimate of 
L using d,nax will be greater than or equal to the number of 
hops required to reach a distance of LR when the differences 
in dn 



Fig. 6. An illustration of i 



(a) 



(b) 



Let A be the area of the intersection of the two radius R 
circles in Fig. 6 and we have from [18] that 

'R-h' 



A = 2\ R^cos 



R 



{R- h)y/2Rh- 



Since A contains N nodes, we have that 

A = N/p. 



(14) 



(15) 



From (14) and (15) we can numerically determine h thus 
giving us 

dmax =R~2h. (16) 
As a result, we need L to satisfy 



R+iL~ l)d,„a. > LR 



which means that 



L = 



R{L - 1) 
R-2h 



1 



(17) 



2) Comparison to Type I Networks: We will compare the 
Type II network simulation results to the Type I analytical 
results. This comparison will allow us to carry over the 
intuition regarding spatial averaging that we have developed 
for the basic cooperative network. However, Type I and Type 
II networks differ primarily in that Type I networks assume 
that all nodes will observe N neighbors from the previous hop 
while any node in a Type II network will only see at least N 
nodes. Thus, if we want to compare Type I and Type II plots, 
we need to establish some meaningful choices of the number 
of cooperating nodes for use with expressions (12) and (13). 

Looking at (a) of Fig. 7, we see that if a node ki in the 
region occupied by the Rk nodes is at the circular boundary 
farthest from node 1 (outer circular boundary), then it will 
likely hear only N nodes from Rk-i- That is, there are TV = 
Alp nodes in area Ai. Recall that N is the minimum number 
of Rk-i nodes any node ki will hear However, looking at 
(b) in Fig. 7, a node ki at the circular boundary closest to 
node 1 (inner circular boundary) in the Rk region will hear 
many more nodes. In fact, a node ki at the boundary between 
Rk~i and Rk will hear the largest average number of nodes 
Nmaxik) = A2p. Since N and Nmax{k) is the range of the 
number of cooperating nodes seen by a node in Rk, it would 
make sense to plot Type I expressions (12) and (13) using 
these two values. However, Nmax{k) varies with k. In Fig. 8 



Fig. 7. (a) Node ki at the outer circular boundary of the Rk set of 
nodes, (b) Node ki at the inner circular boundary of the Rk set. 



we illustrate the regions occupied by the Rk nodes for fc = 1, 
fc > 1, and fc >> 1 overlayed on top of each other and in each 
situation, we see that the set of nodes in Rk seen by a node 
at the boundary between the Rk nodes and the Rk+i nodes is 
different for changing values of fc. However, it is clear that the 
area of intersection always falls inside a semicircle of radius R. 
As a result, we will approximate N,nax = niaxfc:fc>2 Nmax{k), 
by upper bounding the maximum area of intersection with the 
area of the semicircle. This means that 

N,nax~p—- (18) 

Thus, in comparing Type II and Type I results, we will use N 
and Njnax in (18) with both (12) and (13) 



Rjj nodes, k>l 




node (k+l)i 



nodes, k»l 



Fig. 8. The regions occupied by the Rk nodes for fc = 1, A; > 1, and 
fc > > 1 overlayed on top of each other. The region of nodes seen by a node 
at the inner circular boundary of Rk+i changes with fc. 

Using N with (12) and (13) will provide a curve that tends 
to be higher than the Type II simulated curves for two main 
reasons. First, since N is the minimum number of nodes in 
Rk-i that a node ki in Rk will hear and we know that a 
larger number of cooperating nodes will result in decreased 
estimation variance, the variance values computed using N 
will tend to be higher. Second, even if a node ki in Rk hears N 
nodes from Rk-i, each of those N nodes did not necessarily 
only hear N nodes from Rk-2- Thus, the skew and offset 
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estimates made by each of those N nodes in Rk-i whose 
transmissions are being heard by node ki may have a variance 
that is less than predicted by (12) and (13) using N. The 
improved skew and offset estimates made by the nodes in 
Rk-i will thus lead to a lower estimation variance for node 
ki even though node ki hears only N from Rk-i- 

Using Nmax with (12) and (13) will provide a curve that 
tends to be lower than the Type II simulated curves for two 
similar reasons. First, since Nmax is the average number of 
nodes heard by a node at the inner circular boundary of Rk, 
k > 2, and all other nodes in R^ will on average hear fewer 
nodes, a Type I curve using Nmax will tend to yield lower 
values. Second, not all nodes in Rk-i make their estimates 
using a signal cooperatively generated by Nmax nodes. In fact, 
most nodes in Rk-i observe fewer than Nmax nodes. As a 
result, the lower quality estimates made by some of the Rk-i 
nodes will cause the estimation variance of the Rk nodes that 
hear Nmax from Rk^i to be greater than predicted by (12) 
and (13) using Nmax- 

3) Synchronization Performance and Node Density: The 
third issue we want to address in analyzing a Type II network 
deployment is how to decrease synchronization error when we 
know from Section V-B.l that the number of hops L required 
to communicate timing information from node 1 to the edge 
of the network a distance LR away is determined by N. 
Given a fixed R, we can start with some N and p. Using 
(14), (15), and (16), we can determine the value of dmax and, 
hence, from (17) the number of hops L required to send timing 
information from node 1 to the edge of the network. In order 
to decrease synchronization error at a distance LR from node 
1, we need to increase N. However, only increasing N will 
decrease dmax and increase L. Therefore, we need to increase 
both N and p. From (14) and (15), we see that if N/p is kept 
constant, then h wiU be constant. If h is constant, then so is 
dmax- As a result, by increasing node density, we can increase 
the minimum number of cooperating nodes TV and therefore 
decrease synchronization error. 

C. Simulation Results 

In the following simulation results, we have assumed that 
all nodes in the network have no clock skew, i.e. a; = 1 
for all i. From Section IV-C we know that general ai values 
result in variance curves that follow the trends established by 
curves generated using = 1. As a result, using ai = 1 
for all i allows us to study the benefits of spatial averaging 
without considering effects that are dependent on the particular 
network realization. 

1) Comparison to Type I Results: To being the study of 
cooperative time synchronization in general networks, we 
deploy a network for Simulation 1 with the parameters in 
Table II. The simulation results are displayed in Fig. 9. In each 
run, a new network of nodes was uniformly deployed over a 
circular area of radius LR = 5 and the MATLAB simulator 
implemented the cooperative time synchronization protocol. 
Besides plotting the Type I comparison curves described in 
Section V-B.2, we also plot the sample variance of the best 
performing node and the worse performing node. In each 



TABLE II 
Simulation 1 Parameters 



p 


19.10 


N 


4 


R 


1 


L 


5 


d 


2 


m 


4 


a 


0.01 


Number of Runs 


5000 



run, the node in Rk that sees the fewest number of nodes 
from Rk-i is considered the worse performing node while 
the node in Rk that sees the largest number of nodes from 
Rk-i is the best performing node. For the /th run, the fewest 
number of nodes seen by a node in Rk is denoted 
while the largest number of nodes seen by a node in Rk is 
denoted Xmax{k). The skew and offset estimate of the best 
and worst performing node at each hop is recorded and the 
sample variance over the 5000 runs is plotted. 

The top figure in Fig. 9 illustrates the sample skew variance 
curves of the worst and best synchronized node along with 
the Type I curves for comparison. The bottom figure in Fig. 9 
illustrates the clock offset estimate sample variance. Note that 
using equation (17) and the parameters in Table II, we find 
that L = 7- From the simulations, we also see that 7 hops are 
required to traverse the network. In fact, only 7.32% of the 
networks required more than 7 hops to reach all nodes in the 
network. 

As predicted in Section V-B.2, we clearly see in Fig. 9 
that the worst case variance and the best case variance are 
sandwiched between the Type I comparison curves. Also, as 
expected, the skew and offset variances do not closely follow 
the upper and lower Type I comparison curves. The worst 
case skew and offset variance follow the upper comparison 
cruve for the first 2 hops and then begin do deviate from 
the curve. As mentioned in Section V-B.2, this is because 
the nodes contributing to the worst performing node may 
have received signals from more than N nodes. Similarly, 
the best case skew and offset variance follow the lower 
comparison curve for the first 2 hops before deviating. This 
is because many of the nodes contributing signals to the 
best performing node made their estimates using a signal 
cooperatively generated by less than Nmax nodes. Also of 
interest is the steep decrease in the worst case skew and 
offset variance at hop fc = 7. This is due to the fact that on 
average, the distance from the outer circular boundary of the 
i?6 region to the network boundary is much less than dmax- As 
a result, the Rj region is smaller and Xmin (7) will be larger 
than N. Table III shows the - ^ ^™L(fc) 

and Xmax{k\= 5MJ X]f=i° ^™«^('^) values and we see that 
Xmi7i{Q) = N = 4, but Xmin{7) IS nearly twice Xmzn{6)- 

We also note that Xmax{k) increases from 27.56 for fc = 2 
to 35.32 for k = 7. Using (18), however, we find that Nmax ~ 
30. The reason Xmax{k) increases with each hop and does 
not equal Nmax is because Xmax{k) is a different statistic. 
Nmax approximates the average number of Rk^i nodes seen 
by a node ki at the inner circular boundary of Rk- However, 
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Clock Skew Estimate 



o Sample Variance (Worst) 
Upper Comparison Curve 

e Sample Variance (Best) 
Lower Comparison Curve 




Fig. 9. Simulation 1. Top: Sample variance for the skew estimate 
of a Type II network along with Type I comparison curves. Bottom: 
Sample variance of the offset estimate along with Type I comparison 
curves. 

TABLE III 
^min{k) and Xmax{k) for Fig. 9 



k 






! 


1 


1 


2 


4.00 


27.56 


3 


4.00 


29.36 


4 


4.00 


31.86 


5 


4.00 


33.50 


6 


4.00 


34.60 


7 


7.77 


35.32 



of nodes seen by a node ki. Thus, the maximum number of 
nodes would tend to be larger. Note that, not considering the 
effects at the network boundary, X,nin{k) = N because the 
definition of the protocol specifies the minimum to be N and 
there is little randomness in determining Xmin{k)- 
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Clock Skew Estimate 



o Sample Variance (Worst) 
Upper Comparison Curve 

e Sample Variance (Best) 
Lower Comparison Curve 




Fig. 10. Simulation lb. The Type I comparison curves and the Type 
II sample skew and offset variance curves are lower as compared 
to Fig. 9 when TV and p are increased. More cooperation yields 
improved synchronization performance. 



TABLE IV 
Simulation lb Parameters 



Xmax {k) is the largest number of nodes seen by any node ki 
in Rk for the hh network realization. Therefore, Xmaxik) is 
actually an ordered statistic since it takes the largest number 
of nodes seen by a node at hop k. Xmax{k) is thus the mean 
of the ordered statistic. Therefore, we would not expect Nmax 
and Xmax{k) to be the same. Also, X„iax{k) increases with 
k since as the circumference of the circular ring occupied by 
Rk increases, there are more nodes at the boundary between 
Rk and Rk-i- Since there are more nodes at the boundary, 
there are also more opportunities to find the largest number 



p 


23.87 


N 


6 


R 


1 


L 


5 


d 


2 


m 


4 


a 


0.01 


Number of Runs 


5000 



2) Synchronization Performance and Node Density: Next, 
we want to improve synchronization performance by increas- 
ing node density. Starting with the parameters for Simulation 
1, we increase the minimum number of cooperating nodes 
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to iV = 6 while keeping N /p = 0.25 constant. Therefore, 
for Simulation lb (Table IV), p = 23.87 and we plot the 
simulation results in Fig. 10. Comparing Fig. 9 and Fig. 10, 
it is clear that Fig. 10 yields improved skew and offset 
variances, thus showing that increased node density and larger 
N values indeed improve synchronization performance in 
Type II networks. In Table V we show X,„i„(fc) and Xmax{k) 
for Fig. 10. Note that there is only a slight decrease in the worst 
case skew and offset variance curves at hop k = 7 since in 
this simulation, we have that X„ii„(7) = 6.57 is only slightly 
larger than X,nin{Q) = N = 6. 



Clock Skew Estimate of Test Node 



<> Sample Variance 
-A- Upper Comparison Curve 
-y- Lower Comparison Curve 



TABLE V 

T,(k) and Xmaxi^) for Fig. 10 



k 






1 


1 


1 


2 


6.00 


34.01 


3 


6.00 


34.64 


4 


6.00 


37.64 


5 


6.00 


39.50 


6 


6.00 


40.80 


7 


6.57 


41.70 



Another very effective way to visualize how increasing p 
and N can decrease skew and offset variance is to choose 
one test node in the network and consider how its skew and 
offset variance decreases as the network density and number of 
cooperating nodes are increased. In Simulation 2 (Table VI), 
we placed a test node at distance LR = 2.2 from node 1 and 
simulated its skew and offset variance as we increased p and 
N. N took on values ranging from 1 to 10 and we adjusted p 
accordingly to keep N / p = 0.15 fixed. The results are plotted 
in Fig. 11 and we clearly see that as N increases along with 
p, the skew and offset variance of this test node decreases. 
Also, from Section V-B.3, we know that since we keep N/p 
constant, the number of hops required to reach the test node 
stays the same as we increase N. Therefore, since the test node 
is at Z = 3 for every value of N, we have also plotted the 
upper and lower Type I comparison curves for the skew and 
offset variance at hop k ~ 3 to illustrate how the comparison 
curves change in relation to the simulated variance curves. In 
Fig. 11, the simulated skew and offset variance curves of the 
test node fall between the Type I upper and lower comparison 
curves. 

TABLE VI 
Simulation 2 Parameters 



N/p 


0.15 


N 


[1 2 4 6 8 10] 


R 


1 


L 


2.2 


d 


1 


m 


2 


a 


0.01 


Number of Runs 


5000 



It is clear that by keeping the ratio iV/p constant while 
increasing N and p allows us to reduce the synchronization 
error at each hop while keeping the number of hops required to 
synchronize the network, Z, constant. The variance of the skew 




Fig. 11. Simulation 2. Variance of the skew and offset estimates of 
the test node fall between the Type I comparison curves and decrease 
with increasing N and p. 



and offset estimates is decreased by increasing the minimum 
number of cooperating nodes. 

Furthermore, from the simulations in this section, we find 
that the upper and lower Type I comparison curves provide a 
good reference to the performance of Type II networks. We 
have established that the best and worst case variance values 
for the Type II skew and offset estimates fall between the 
upper and lower Type 1 comparison curves. As the density of 
the network and iV are both increased, the comparison curves 
will shift downwards and become closer together Thus, we 
would expect the variance of the Type II network estimates to 
change similarly with increasing N and p. 

VI. Conclusion 

In this paper we have proposed one technique that uses 
spatial averaging in dense networks as a means to improving 
global time synchronization. Spatial averaging is used to 
improve the timing data points that are used to estimate 
clock skew and clock offset. By decreasing the error in the 
timing data points, improved clock skew and clock offset 
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estimates can be made. Our analysis of the technique in a basic 
cooperative network revealed that the error variance in both 
the clock skew and clock offset estimates can be significantly 
decreased as the number of cooperating nodes increases. Simu- 
lation results also show that synchronization over large, multi- 
hop networks can be improved by increasing node density. 
Further analysis and a comparison between cooperative and 
non-cooperative techniques can be found in [19]. 

This scalable protocol provides an alternate way to combat 
the scalability problem. It allows us to simply increase the 
number of nodes in the network to obtain improved syn- 
chronization performance. The new trade-off between network 
density and synchronization performance provided by spatial 
averaging will allow for added flexibility in designing future 
networks. 

It is important to note that the concept of spatial averaging is 
very general and our proposed cooperative technique is but one 
manner in which to take advantage of it. Our protocol shows 
that techniques using spatial averaging can be designed. Even 
though the proposed protocol has certain limitations, such as 
requiring access to the physical layer, it allows us to success- 
fully illustrate the performance improvement achievable using 
spatial averaging. Future work will focus on other approaches 
to spatial averaging. For example, it would be desirable to 
develop a cooperative technique using spatial averaging that 
achieves performance gains while needing only access to the 
data link or network layer 

Appendix 

Proof of Theorem 1 Node 1 begins the synchronization 
processes by transmitting a sequence of pulses at times tq + W, 
for / = 0, . . . , 771 — 1. For simplicity, assume that tq and d are 
integer values. Note that since node 1 transmits these pulses 
in its own time scale ci (the reference time), the pulses will 
occur at integer values of t. Using the clock model in (1), any 
node li, z = 1, . . . , iV, in the Ri set of nodes will get a vector 
of observations Yi^, where Yii[l] = aii(To — Ai,) + 
and the (/ + l)th element of Yi^ is Yii[l + 1] = aii(ro — 
Aii) + Idaii + ^Pii f+i. This can also be written as 



where 



Yi, = 













(19) 



aii(To - Aii) 
an 



with H as in (4) and = [Wus, ■ ■ ■ ,Wu^rnV . Since 

is an independent Gaussian random variable for each 
I, Wii ^ A/'(0, Sii) with Eli = c^Im- As mentioned, this set 
of observations is for any node li in the set of Ri nodes. 

Since we have N Ri nodes, we can write the vector of 
observations made by all Ri nodes as 

(20) 







Yi = 


H01 + 


Wi 


( 


where 


" Yii " 




" ^11 " 




" Wn ' 


Yi = 


. YiTv . 


, Oi = 




, Wi = 


. Wi^ . 



and H is as in (6). Note that Wi - A/'CO, cr^I^,,,). This 
way we have Yi as the vector of observations made by all 
Ri nodes and we can make a UMVU (uniformly minimum 
variance unbiased) estimate of 6i by taking 



where 



(H^H)-iH^Yi ~A/-(/ii,Ei) 
^1, Ei = a2(H^H)-\ 



Ml 



It is easy to see that 
(H^H)-i 



(H^H) 



(H^H)- 



and 



(H^H)-i 



2(2tTt-l) -6 

m{m-{-l) dm(m-l-l) 

~6 12 

dm{m+l) l)m(m+l) 



This establishes the initial conditions for the theorem. 9i is a 
2N X 1 column vector where the subvector made up of the 
(2(i - 1) + l)th and (2(i - 1) + 2)th elements, i^l,...,N, 
is dii = (H^H)~^H^Yii. Therefore, any node li's skew 
estimate (2) and offset estimate (3) can be found from 9i as 



"li 



'2(i-l)+2'^l 



and 



(21) 



(22) 



where e/ is the column vector of all zeros except for a one in 
the kh position. 

Each node li can now make an estimate of the next appro- 
priate integer value of t, in this case t = tq + md, by making 
a minimum variance unbiased estimate of 6'ii i + mdOii^2 = 
aii{To — Aii) + mdau- This can be done with the estimator 



mdOii 



where Co = [1 md\. This will then be node li's estimate of 
the next appropriate integer value of t in its own time scale 

Cli. 

From (5), every node li will then transmit a sequence of 
TO pulses occurring, in the time scale of cii, at Xii{l) ~ 
Tii + ld6ii,2, for I ~ 0, . . . , TO — 1. Using the clock model 
(1), we find that in the time scale of ci these pules occur at 



(fij + ld9u^2)ci 



IdOu 



an 



an 



Id- 



'li,2 



+ An 

*liJ+l 



an an 

Any node 2j in the R2 set of nodes that can hear node li will 
thus get a sequence of pulses 



Y 



Tn 

an 



An + Id- 



ni.2 



*i 



an 

+ *2j,/+l, 



an 
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where Z = 0, . . . , m — 1. 

In this Type I network deployment every node 2j hears the 
same set of N nodes and takes the sample mean of each cluster 
of pulses for its observation, we can express the actual vector 
of observations made by node 2j as 



where 



N 



an 



Id- 



'li,2 

an 



au 



2i,/+i7 



where Z = 0, . . . , m— 1. Note that since these pulse arrivals are 
clustered, we assume that for a given cluster, each pulse arrival 
is corrupted by the same jitter. Thus, receiver side jitter ^'2j.i+i 
is an independent sample for every I, but takes the same value 
for each i. This models the fact that clock errors occurring in 
a small time window are highly correlated while errors farther 
apart in time are independent. We can rewrite this simply as 
Y2j[^ + 1] = a2j((ri+/dai-*i,;+i)-A2j) + 'I'2i,i+i, where 



TV „ 

2—1 



Ah 



- A i \ - 

2 — 1 



li,2 



N 



*1,( + 1 = Tr 



an 



Since every node 2j will see the same N, this means that 
every node 2j will have the same ti and ai. Therefore, ri 
and ai are now fixed, and it can be easily found that 
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Node 2j's vector of observations can also be written 


linear form similar to (19), Y2j 


= H02, - 


fW2j, 


where 
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The vector of observations made by all R2 nodes can be 
written in a manner similar to (20), 



Y2 



H6I2 + Wo 



" Y21 " 




O2I 




a2llm 




, 02 = 




, Q2 = 




Y2^ 




^2N 




a2N'^m 



W2 



W21 
W22 



w 



2N 



Q2 



*1,1 
. ^l,m 



*21,1 
*21.m 
*2W,1 
^2N.m 



This means that W2 ~ A/'(0, S-^^), where 
The i?2 nodes will estimate O2 as 

I2 = (H^H)-iH^Y2 (23) 
~ AA(0-2,(H^H)-iH^I]w,((H^H)-iH^f). 

However, for analysis, this does not give us the complete 
distribution of O2 since O2 is a function of 9i. Therefore, we 
first consider how 62 j is a function of 61. We find that 

a2j{Ti - K2j) 

a2ja\ 

«2,(^Eillt7+^l''-^2,) 

0-2] 2^i=i an 
, I s-^N Bii i+md9ii 2 I A A \ 



^2i 



"2j ^ Z^i=l 



a2j 
iV 



TV 

E 

2 = 1 





1 dm ' 






( 


an aii 



L an J 




Oli.2 



Ah 




"2j 



A2, 





Using (24) we have 



92 — A.201 + B2 



(24) 



(25) 



where 



D2 



1 dm 
Qii ail 

^ 

ail 















dm 



B2 = D2 



" An " 




a2iA2i 














'^2wA27V 
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for 



N 



a2i 
a2i 



a 



2N 

a. 







2N 



a2i 
a2i 



2JV 

a. 







2Ar 



Using (23) and (25), the distribution of 02 can now be found. 

122 = E{§2) 

= E{E{e2\k)) 

= Eih) 

= £;(A2^i+B2) 

a2i(To + md - A21) 

^21 



(26) 



"2W 

Using the decomposition 

Cov(^2) =i;(Cov(l2|^i)) +Cov(£;(^2|^i)), 
we have from (23) and (25) 

= E{Cow{e2\ei)) 

= (H^H)-iH^Sw,((H^H)-iH^)^ 



Cov(£;(6i2|6ii)) = Cov(6i2) 



giving us 



S2 = Cov(^2) = 



A2SiAf. 



(27) 



Thus, the distribution of 62 is 

4 -A/'(M2,S2). 

^2 is again a 27V x 1 column vector where the subvector made 
up of the (2(i - 1) + l)th and (2(i - 1) + 2)th elements, i ^ 
l,...,iV, is §21 = (H^H)-iH^Y2,. Therefore, as in (21) 
and (22), any node 2i's skew estimate (2) and offset estimate 
(3) can be found from 62 as 



and 



a2i - e2(j_i)+2P2 



A2i = ej(,_i)+i6'2 - (to + dm). 



(28) 



(29) 



Each node 2i will now be able to transmit a sequence of 
TO pulses occurring, in the time scale of C2i, at X2i{l) = 
T2i + ld92i,2, for / = 0, . . . ,TO-1, where f2j = 02is+md02i,2- 
Repeating the same process we carried out for the observations 
of any node 2j with any node 3j, we can find 63 ~ J^ip-s, S3). 
In fact, continuing this procedure, we can find the distribution 
of 6k for the Rk nodes as 



where similar to (26) we have /ifc = E{0k) which is found in 
(7) and similar to (27) we have 



tk = Cow{h 



k I 



which is found in (8). S„ifc> A^, B^., and S-^^ are as in the 
theorem statement. As in (28) and (29), any node fcz's skew 
estimate (2) and offset estimate (3) can be found from 9k as 



aki 



°2(i~l)+2"k 



and 



Afcj = ej(i_i)+i6'/c - (to + dm{k - 1)). 
This concludes the proof of Theorem 1 . A 
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